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This paper presents both a numerical method for general relativity and an application of that 
method. The method involves the use of harmonic coordinates in a 3+1 code to evolve the Einstein 
equations with scalar field matter. In such coordinates, the terms in Einstein's equations with the 
highest number of derivatives take a form similar to that of the wave equation. The application is 
an exploration of the generic approach to the singularity for this type of matter. The preliminary 
results indicate that the dynamics as one approaches the singularity is locally the dynamics of the 
Kasner spacetimes. 
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I. INTRODUCTION 

It has long been known that harmonic coordinates are 
useful for mathematical relativity. In particular, these 
coordinates were used jjj to prove local existence of so- 
lutions of the vacuum Einstein equation. This usefulness 
of harmonic coordinates stems from their putting Ein- 
stein's equation into a form that is similar to the curved 
spacetime wave equation. Since many numerical tech- 
niques work well on the wave equation, one might expect 
that harmonic coordinates would be used extensively in 
numerical relativity, and it is somewhat surprising that 
they are not. (However see Q for some recent mathe- 
matical and numerical work on the linearized case. In 
addition, harmonic time slices have been advocated and 
used in numerical relativity ||-0). This is perhaps due 
to the following drawback of harmonic coordinates: these 
coordinates are solutions of the wave equation, and such 
solutions need not have a timelikc gradient at all points 
of spacetime, even if they start out with a timelike gradi- 
ent on an initial data surface. What this means is that in 
harmonic coordinates, the time coordinate will in general 
not remain timelike and this is likely to cause numerical 
problems. |||| As we will see later, there is a way around 
this difficulty. 

One area of numerical relativity where harmonic coor- 
dinates have been used (though somewhat unintention- 
ally) is in the study of the approach to the singularity, in 
particular in the Gowdy spacetimes. fic|| Numerical sim- 
ulation of approach to a spacetime singularity presents 
problems of its own. Since it is expected that various 
quantities become infinite at the singularity, the numeri- 
cal simulation will generally stop after a finite coordinate 
time, the time when a surface of constant time first en- 
counters the singularity. Since this first encounter gen- 
erally occurs at one spatial point, the information about 
the behavior of the singularity at other spatial points 
will be unavailable from the numerical simulation. The 
solution to this difficulty is to choose a time coordinate 
that tends to infinity as the singularity is approached. 
In this way, the simulation is not forced to end at finite 



coordinate time, the whole spacetime up to the singu- 
larity is covered and the behavior of the metric as the 
singularity is approached simply becomes asymptotic be- 
havior in the limit of large time coordinate. For Gowdy 
spacetimes there is a natural choice of such a time co- 
ordinate: these spacetimes are foliated by T 2 s invariant 
under the symmetry group. The area of the symmetry 
T 2 s goes to zero as the singularity is approached. There- 
fore minus the logarithm of this area is a natural time 
coordinate that goes to infinity as the singularity is ap- 
proached. While this method works well for the Gowdy 
spacetimes, since it depends crucially on the symmetry 
of the Gowdy spacetimes, the method does not seem to 
generalize to the case of the generic singularity with no 
symmetries. A natural generalization comes when one 
notices that this time coordinate (minus the logarithm 
of the area of the symmetry T 2 s) is also harmonic. Since 
in some sense one expects the wave equation to become 
singular as the spacetime singularity is approached, one 
might also expect a solution of the wave equation to blow 
up as the singularity is approached and thus one might 
want to use such a solution as the coordinate time. 

What is expected to be the generic behavior of a space- 
time as the singularity is approached? Based on studies 
of spacetimes with T 2 symmetry |Io|-|l2[| and spacetimes 
with £7(1) symmetry [Q, the expected answer [|l4| is 
the following: the singularity is expected to be spacelike 
and as it is approached each spatial point is expected to 
"decouple" from the others and undergo a dynamics cor- 
responding to that of a homogeneous spacetime (though 
a different homogeneous spacetime at each spatial point). 
Which types of homogeneous spacetime is this dynamics 
expected to correspond to? For vacuum spacetimes, it 
is thought that the dynamics will be oscillatory, possibly 
corresponding to the Mixmaster spacetime as conjectured 
in p| . For many other types of matter it is expected that 
as the singularity is approached the matter terms in the 
Einstein equations become negligible and the dynamics 
approaches that of a vacuum spacetime. 

However for so called "stiff matter" (i.e. a scalar field 
or a perfect fluid with equation of state P = p) it is ex- 
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pected that the dynamics will not be oscillatory and will 
correspond to that of a Kasner spacetime. This expecta- 
tion is greatly bolstered by a theorem due to Andersson 
and Rendall jl6| which shows local existence in a neigh- 
borhood of the singularity of solutions to the Einstein- 
scalar equations, with the expected asymptotic behavior 
and with enough degrees of freedom to be the generic 
solutions. 

Since the approach to the singularity is expected to be 
simpler for stiff matter than for vacuum, the stiff matter 
case should be easiser to treat numerically. Therefore, in 
this paper we confine ourselves to a numerical study of 
the approach to the generic singularity in the Einstein- 
scalar system. Section 2 presents the equations and nu- 
merical methods used. The results are given in section 3. 
Section 4 contains a discussion of the results and of other 
possible applications of the harmonic coordinate method. 

II. EQUATIONS AND NUMERICAL METHODS 

The equations that we wish to evolve numerically are 
the Einstein-scalar equations 

R al3 = 8ttV a <t>V (f> (1) 

Here, we use the conventions of |l7|] including units where 
c = G = 1. As a consequence of Eq. ([!]) and the Bianchi 
identities, the scalar field must satisfy the wave equation 

V Q V> = (2) 

The Ricci tensor is given in terms of the Christoffel sym- 
bols by 

R a p = d-r^ap — da^p + r^r^ 7 — r^r^ (3) 

while the Christoffel symbols are given in terms of the 
metric by 

T lp = h lS {d a 9f3S + dpg aS - d s g a p) (4) 

Harmonic coordinates are solutions of the wave equa- 
tion. As a generalization of harmonic coordinates, con- 
sider coordinates that satisfy the wave equation with 
source 

V a V a x" = H» (5) 

where H 11 are specified from the beginning. Then using 
Eqs. (||(|,||) we find that the Ricci tensor is given by 

R a p = — \9 1 ° d 1 d (y g a p + ^Cp^C^a 
+ \C a ^C^ - T2 a T^ - d (a H p) + Tl p H, (6) 

where C a ^ v = d a g^ v . Note that the second derivative 
terms appear only in the wave operator. Therefore, one 
might expect that Einstein's equations in this form be- 
have similarly to the wave equation and that numerical 
methods that work well on the wave equation might work 



well on Einstein's equations in this form. The reason for 
considering nonzero source terms H» in Eq. (|) is that 
these terms allow us to change the behavior of the time 
coordinate and thus may allow us to eliminate (or at least 
postpone) the behavior where the time coordinate ceases 
to be timelike. In the simulations done in this paper, no 
source terms were needed. I know of no systematic way 
to find appropriate source terms and expect that in the 
cases where they are needed, one will have to resort to a 
trial and error method to find appropriate H^. Note that 
the use of spatial harmonic coordinates could also lead to 
coordinate problems. This would occur if the gradients 
of the four coordinates fail to be linearly independent. If 
this sort of problem occurs, one might expect to be able 
to postpone or eliminate it by using appropriate source 
terms. 

The numerical method used requires equations that are 
first order in time. To put the equations in such a form, 
we define quantities P a/ 3 and P<j> given by 

Pap = d t g a f3 (7) 

P+ = d t ct> (8) 

Then the Einstein-scalar equation becomes 

- g W d t P a p = 2g ok d k P aP + g ik did k g a(} 
+16Trd a cf ) d cj ) + 2d (a H f3) - 2Tl /3 H 1 
-CcTC^p - CfTC^ + 2TZ a T»p (9) 

The wave equation for (f> becomes 

- g oa d t P^ = 2g ok d k P^ + g^didkcj* - g aS3 V\d^ (10) 

The full set of equations that are evolved in the computer 
code are Eqs. ([t|-|To|) . 

We now turn to the numerical methods used to evolve 
these equations. Spatial derivatives are approximated by 
centered differences. The variables are evolved in time us- 
ing a three step iterated Crank-Nicholson (ICN) method. 
p"8| p0| This works as follows: evolution equations of the 
form d t S = W(S) for some set of variables S are approx- 
imated as 

S n+1 = S n + ^[W(S n ) + W(S n+1 )] (11) 

where S n is the value of S at time step n and At is the 
time step. Then, using S n as an initial guess for S n+1 , 
Eq. (|Ti"| ) is iterated three times. 

The spacetimes we consider have topology T 3 x R. 
Each spatial slice has topology T 3 . In terms of the spa- 
tial coordinates, this means that < x < 2ix with and 
2ir identified (and correspondingly for y and z). This 
topology is implemented numerically as follows: a spa- 
tial coordinate x has TV grid points. The variables on 
points from 2 to N — 1 are evolved using the evolution 
equations. The variables on point 1 are set to the values 
at point N — 1 while at point N they are set to the values 
at point 2. 
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III. RESULTS 

All runs were done in double precision on Compaq 
XP1000 workstations and on the NCSA Origin 2000. 
The time step was At = Ax/2. While the source terms 
should be helpful in keeping the coordinates well be- 
haved, we did not need them for the cases studied here 
and all runs have = 0. 
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FIG. 1. comparison of the quantity P as found using a 
Gowdy code and the 3+1 harmonic code 



code, 3 gridpoints were used in the x direction, 3 grid- 
points in the y direction and 500 gridpoints in the z 
direction. For the comparison the initial data used is 
P = 0, d r P = 5cosz, Q = cosz, d T Q = 0. These data 
are evolved until t = ir and the results for the compar- 
ison are given in Fig. 1. Here, the solid line represents 
the 1+1 evolution of the Gowdy equations, while the dots 
represent the full 3+1 evolution using the harmonic code. 
There is clearly agreement between the two. 
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FIG. 2. convergence test involving the constraint 



Before exploring the generic singularity, we would like 
to test the code. One way of doing this is to run the 
code on cases with symmetry where results are already 
known. Such a case is the Gowdy spacetimes. These 
have the form 10 



ds 2 = e^' 2 [~e- 2T dr 2 + dz 



+e- T [e p dx 2 + 2e p Qdxdy + (e p Q 2 + e - p ) dy 2 ] (12) 

Here, P, Q and A are functions of r and z. It follows from 
Eq. ( |l2] ) that the coordinates (r, x, y, z) are all solutions 
of the wave equation. Therefore, they are harmonic coor- 
dinates. The vacuum Einstein equations for the Gowdy 
spacetimes are |ll]] 



d T d T P - e- 2r d z d z P - e 



■if 



(cVQ) 2 - e~ 2T 



{dzQY 



= 



(13) 



Before presenting the results of another code test, we 
turn to the question of finding initial data without sym- 
metries. On an initial data surface, the intrinsic metric 
hij and the extrinsic curvature Kij must satisfy the con- 
straint equations 



DilCj - DjK = Sii^Dji 



Mr- 



K 2 - K ll K %j = 



D l (f)Di 



(15) 



(16) 



(Note that we use the convention of |l7j for the sign of 
which is opposite to the convention usually used in 
numerical relativity). Here, an overdot denotes deriva- 
tive along the normal to the surface, and Di and 
are respectively the covariant derivative and scalar cur- 
vature of hij. Given a solution of Eqs. (|l5l ) and (|l(|), 
we produce initial data for evolution in harmonic coor- 
dinates as follows: For spatial directions i and j we have 



hij, g%o — 9oi — 0, goo 



-1, P, A 



2K li . The 



d T d T Q - e 2r d z d z Q + 2 (d T Pd T Q - e 2r d z Pd z Q) = remaining components of P af3 are solved for using Eq. 



(14) 

plus constraint equations that determine A once P and 
Q are known. 

To test the 3+1 code, we evolve Eqs. (|l3|) and |L4|) 
using a 1+1 code. Then we evolve the same initial 
data with the 3+1 harmonic code and compare the re- 
sults. To do the comparison, note that P = r + In g xx 
so that it is straightforward to compare the values of 
P produced by the two codes. For these simulations 
500 gridpoints were used in the 1+1 code. In the 3+1 



We want to find a solution of Eqs. ( |l5| ) and ( IHf ) that is 
simple but has no symmetries and has some free parame- 
ters. We choose = and hij equal to the flat Euclidean 



metric in the usual coordinates, 
ture we choose 



K 



For the extrinsic curva- 



(bi + a 2 cos y + a 3 cos z) /2 
K yy = (62 + cl\ cos x — 03 cos z) /2 
K zz = (b 3 — ai cosx — a2 cosy) /2 
K xy = K yx = (ci cos z) /2 
K xz = K zx = (c 2 cos y) /2 
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K 7 



(C3 cosx) /2 



(17) 



Here, the quantities a 2 ;, bi and q are constants that are 
free parameters. It is straightforward to show that the 
extrinsic curvature of Eq. (|17| ) with our choice of initial 
hij and (f> satisfies Eq. (|15|). Equation ( [l6|) then be- 
comes an agebraic equation for which can be solved 
provided that the left hand side of the equation is pos- 
itive. Note that if all the ai and Cj are zero and all the 
bi are equal, then the initial data evolve to the spatially 
flat Robertson- Walker spacetime with scalar field matter. 
Thus, this family of data can be thought of as Robertson- 
Walker with large gravitational and scalar waves. 

We now consider a convergence test involving a con- 
straint that comes from the use of harmonic coordinates. 
Define the quantities by 



y 1 a8 



(18) 



(this is the appropriate form of the constraint for the 
case where the source term vanishes. For the general 
case, the constraint would be given by = g a ^T^ + 
Hi 1 ). Then from Eq. (§) it follows that C = 0. Since 
we are solving the evolution equations by approximating 
them by finite difference equations, the quantities C as 
evaluated by the computer code will not be zero because 
of errors due to the finite grid spacing Ax. Define the 
quantity C by 



C 



J y/gg a i3C a C p dxdydz 



J ^/gdxdydz 



1/2 



(19) 



This quantity is a type of measure of the average size 
of the constraint. Figure 2 shows a plot of C vs. 
time. The parameters used are = (0.1,0.1,0.2), bi = 
(0,-0.5,-0.5) and c$ = (0,0,0). Here the curve corre- 
sponds to a run done with 20 gridpoints in each spatial 
direction, while the dots correspond to a run with 38 
gridpoints in each spatial direction (which gives half the 
grid spacing) and with C multiplied by 4. The results 
show second order convergence. 

We now consider the approach to the singularity. To 
see what is expected, it is helpful to consider the Kasner 
spacetime in harmonic coordinates. This is given by 



ds 2 



. e (<H+Q2+q 3 )r dT 2 + + ^2^2 + ^3 



(20) 



where the qi are constants This metric is generally a so- 
lution of the Einstein-scalar equations; but for the case 
where (X^ft) 2 = it i s a vacuum spacetime. (In 

the usual treatment of vacuum Kasner spacetimes, one 
defines the quantities pi = 5i/Q3<2m.) and then obtains 
the condition J^Pi = YlPi — Note that in the vac- 
uum case one cannot have all three directions contract- 
ing but that this is possible for the Einstein-scalar case. 
Note also that the metric components are exponential 
functions of time. It also turns out that the scalar field 
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FIG. 3. behavior of metric components and scalar field as 
the singularity is approached at the spatial point (0, 0, 0) 




In (-gjt) 



FIG. 4. behavior of metric components and scalar field as 
the singularity is approached at the spatial point (0, 7r/4, 7r/2) 



is a linear function of time. Thus, if the behavior of a 
generic solution near the singularity is local Kasner, then 
we should expect metric components that are exponen- 
tial functions of time, with the exponent depending on 
space. We should also expect a scalar field that is a linear 
function of time with the slope depending on space. 

Results on the approach to the singularity are shown 
in Figs. 3 and 4. The run was done with 34 grid points in 
each spatial directions. The parameters aj, bi and Ci are 
the same as for the convergence test. Here the scalar field 
and the logarithms of the diagonal metric components 
are plotted as functions of time. Figure 3 corresponds 
to the spatial point (x,y,z) — (0,0,0) while for Fig. 4 
the point is (0, 7r/4, 7r/2). Note that as the singularity is 
approached these quantities all become linear functions 
of time. The differences between Figs. 3 and 4 show 
that there is a spatial dependence of the approach to the 
singularity. Note from Fig. 3 that though the x direc- 
tion is initially expanding, eventually all three directions 



4 



contract. This is what one would expect if the metrics of 
||l6f represent the generic behavior near the singularity 
since these metrics have all three directions contracting 
in a neighborhood of the singularity. 

IV. DISCUSSION 

While this study is somewhat preliminary, it indicates 
that harmonic coordinates can be a useful tool in numer- 
ical relativity. Though, in principle coordinate problems 
could occur, this did not happen in the cases studied 
here, even though they involved very strong fields. Fur- 
thermore, the use of the source terms may cure such 
problems if they arise. 

As for the behavior of generic singularities, the numer- 
ical results indicate that solutions of the form proved in 
reference Jl6[ to exist in a neighborhood of the singular- 
ity also exist globally. Thus such solutions are likely to 
describe the generic singularity in the stiff matter case. 

There are several projects for which the methods of this 
paper could be used. One is to do a more extensive study 
of the singularity in the Einstein-scalar case, with a more 
thorough exploration of the evolution corresponding to 
various values of the parameters in the initial data of the 
previous section. It would also be helpful to evolve for a 
longer time. 

Another project is to remove the scalar field and study 
the approach to the singularity of the generic vacuum 
spacetime. Here, the behavior is expected to be more 
complicated and a treatment will probably require more 
spatial resolution to resolve the expected sharp features, 
as well as longer evolution in time to see the expected 
oscillatory behavior. 

Yet another project is to study the behavior of asymp- 
totically flat spacetimes rather than closed cosmologies. 
Here, the closed cosmologies were studied partly for sim- 
plicity. The periodic boundary conditions arc simple 
to implement and completely consistent with Einstein's 
equations. In contrast, for an outer boundary at a finite 
distance in an asymptotically flat spacetime, one needs 
to put some sort of outgoing wave boundary condition. 
Such conditions are usually not consistent with Einstein's 
equation (it is known how to have a consistent boundary 
condition [|T) but this condition is quite complicated). 
These inconsistent conditions may lead to numerical in- 
stability. Since harmonic coordinates make Einstein's 
equation look like the wave equation, simple outgoing 
wave boundary conditions that work numerically with 
the wave equation might be expected to "work" (at least 
in the sense of not causing numerical instability) for Ein- 
stein's equation. 

Since much work has been done on numerical simula- 
tions of asymptotically flat spacetimes using the standard 
Arnowitt-Deser-Misner (ADM) j22) approach, it is help- 
ful to make comparisons with this approach. In the ADM 
approach, the spacetime metric is written as 



ds 2 = -a 2 dt 2 + h l3 {dx l + I3 l dt){dx j + (3 3 dt) (21) 

Einstein's equations are written as an evolution equation 
for the spatial metric hij (and the extrinsic curvature 
Kij) while the gauge choice results in equations for the 
lapse a and the shift f3 l . This framework is sufficiently 
general to accomodate the use of harmonic coordinates, 
which correspond to the following equations for lapse and 
shift 

d t a - frdia = Ka 2 (22) 

d t (¥ - (3 m d m (3 l = h mn r mn a 2 (23) 

(Note that the sign of K in equation ( p2| ) comes from the 
convention of Jl7[] ) - Here T % mn is the Christoffel symbol 
associated with hij. (Equations (p3) and ( p3] ) hold for 
the case of vanishing source term H^. Similar equations 
hold in the case of nonzero iJ M ). The use of equations 
( E^ ) and ( [23| ) in an ADM code is not precisely equiva- 
lent to the approach of this paper. The reason for this 
is that equations ( p2| ) and ( f23| ) directly solve the con- 
straint 5 = 0, while in our approach, this con- 
straint is used to change the form of the evolution equa- 
tions. Nonetheless, it would be interesting to use equa- 
tions (p2|) and (p3| ) in an ADM code to see how that 
compares with other choices of lapse and shift. In par- 
ticular, one might expect better stability properties and 
compatibility with a simple outgoing wave boundary con- 
dition. (Though note that such improvements can also 
be obtained using the BSSN approach |25|-|2"5f| ) . 

Another desired feature of a numerical code is the abil- 
ity to treat black holes, and this sometimes requires black 
hole excision. While harmonic coordinates are singularity 
avoiding || they are just barely so and come arbitrarily 
close to a singularity. Thus the need for excision in an 
approach that uses harmonic coordinates should be at 
least as great as in the standard approach. Nonetheless, 
one might hope that excision itself would be easier to 
implement using the approach of this paper. This is be- 
cause no elliptic equations are involved and the light cone 
of the wave operator is the same as that of the physical 
metric. 

All these projects are work in progress and prelimi- 
nary results from them are promising. Thus, I expect 
that harmonic coordinates will become a useful tool in 
numerical relativity. 
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